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Abstract 

We study molecular dynamics within populations of diflFusivcly cou- 
pled cells under the assumption of fast diffusive exchange. As a technical 
tool, we propose conditions on boundedness and ultimate boundedness for 
systems with a singular perturbation, which extend the classical asymp- 
totic stability results for singularly perturbed systems. Based on these 
results, we show that with common models of intracellular dynamics, the 
cell population is coordinated in the sense that all cells converge close to 
a common equilibrium point. 

We then study a more specific example of coupled cells which behave 
as bistable switches, where the intracellular dynamics arc such that cells 
may be in one of two equilibrium points. Here, we find that the whole 
population is bistable in the sense that it converges to a population state 
where either all cells are close to the one equilibrium point, or all cells 
are close to the other equilibrium point. Finally, we discuss applications 
of these results for the robustness of cellular decision making in coupled 
populations. 

1 INTRODUCTION 

Coordinated behavior in cell populations is a long standing topic in the analysis 
of biological systems. While so far mainly the synchronization of biological 
oscillators has been in the focus of this research [9, 12, 8], coordinated behavior 
is also of interest for other types of cellular behavior. The study presented in this 
paper is motivated by the analysis of biological switches, which are a common 
phenomenon in biological systems such as gene regulation or signal transduction 
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[5, 4]. In contrast to oscillators, synchronization among biological switches has 
hardly been studied so far. Wang and coauthors [14] present a model of a genetic 
switch in a cell population where one of the molecules making up the switch is 
diffusively exchanged among cells via the extracelhilar medium. However, [14] 
is concerned more with a numerical study of the effects of noise on synchronized 
switching than the underlying mechanism of switch synchronization per se. A 
more theoretically oriented stiidy in this area has been performed by Schmidt 
and coauthors [11]. There, conditions for consensus among directly coupled, 
scalar bistable switches are proven using Lyapunov techniques. While no specific 
biological model is studied in [11], the analysis is motivated from the occurence 
of bistable dynamics in biological systems. 

Here, we apply singular perturbation theory to study populations of diffu- 
sively coiiplcd cells for diffusion rates which are significantly faster than the 
intracellular dynamics. Although there are long-standing results for the conver- 
gence analysis of singularly perturbed system [10, 3], these studies considered 
only systems where the perturbation does not affect the position of the equilib- 
rium point, and then proceeded to show asymptotic or exponential stability of 
this equilibrium point for a sufficiently small perturbation. In the model we are 
going to construct here, the equilibrium point will in fact depend on the per- 
turbation, and the previous results are not applicable. We can however prove 
a new result on boundedness and ultimate boundedness of singularly perturbed 
systems, which in general is applicable to cells coupled via diffusive exchange 
of molecules with a common extracellular medium. As a specific application, 
we consider a population of coupled bistable switches. The theoretical results 
and a simulation study show coordinated behavior and synchronized switching 
among the coupled switches for fast enough diffusive exchange rates. 

The paper is structured as follows. In Section 2, we present conditions on 
boundedness and ultimate boundedness of singularly perturbed systems. In 
Section 3, we construct a generic model for a population of diffusively coupled 
cells. The previous conditions are then applied to show convergence to the 
neighborhood of a common equilibrium point for all individual cells. Finally, in 
Section 4, we consider a population of diffusively coupled bistable switches, and 
show that the population behaves as a single switch with all cells in coordination 
for fast enough diffusive exchange. 

2 BOUNDEDNESS IN SINGULARLY PERTURBED 
SYSTEMS 

2.1 Conditions for ultimate boundedness 

The analysis of the coupled switches with singular perturbation requires an ex- 
tension of singular perturbation theory. Previous singular perturbation results 
either assumed that the equilibrium point does not vary with the small parame- 
ter e [10], or that the dependence on e of the dynamics tends to zero as the state 
approaches the equilibrium point [3], and then proceeded to show exponential 



stability for this equilibrium point. In this study, these assumptions arc not 
satisfied and it will not be possible to show exponential stability of a fixed equi- 
librium point. Instead, we will show convergence of the singularly perturbed 
system to a neighborhood of the equilibrium state. By decreasing the small 
parameter e, this neighborhood can be made arbitrarily small, thus recovering 
the classical result in the limit. 
Consider the system 

x = fix,z) 
ez = g{x,z,e) 

with X G Dx C M." , Dx convex, z G C W", e > 0, and functions / and 
g continuously differcntiable in all of their arguments and such that a unique 
solution x{t) S Dx, z{t) € exists for all t>0. Initial conditions are taken as 
a;(0) = xo and z{0) = zq. 

Our first assumption is that a quasi-steady-state solution for the fast dy- 
namics can be obtained, and that the system (1) has an equilibrium point at 

£ = 0. 

Assumption 1. There exists a continously dijferentiable function h : — >• 
such that 

gix,hix),0) = (2) 
for all X G Dx, and h satisfies the linear growth hound 

ll|^(^)ll<ci (3) 

for a positive constant ci for all x G D^- Furthermore, for e = 0, system (1) 
has an equilibrium point at (0, /i(0)), i.e., 

/(o,Mo)) = o 

5(0,/i(0),0) = 0. 

For the next step, the reduced order systems for £ = are defined: the slow 
system for the variable x and the boundary-layer or fast system for the variable 

y{T) = z{r) - h{xo) (5) 

on the fast time scale t = t/s. The slow system is defined as 

X = f{x,h{x)), (6) 

and the boundary-layer system as 

^=g{xo,y + h{xo),0) (7) 
where Xq is considered as a constant parameter. Also, define 



Af{x, z) = f{x, z) - f{x, h{x)) 



(8) 



and 

Ag{x,z,e) = g{x,z,e) - g{x,z,0) (9) 

We have the fohowing assumptions on the reduced order systems and the 
functions / and g: 

Assumption 2. There are a scalar function Vs{x) and positive constants C2, 
di, 6,2, and D such that 

di\\xf <Vs{x) <d2\\xf (10) 

and 

f{x,h{x))<-C2\\xf (11) 



dx 

for all X e 0.x, a compact subset of D^, with {di||a;||^ < D} c O^. 

Assumption 3. There are a scalar function VB{x,y) and positive constants cs, 
ds, and di, such that 

dzMf <VB{x,y)<d4y\f (12) 

and 

g{x,y + h{x),Q)<-c^\\yf (13) 



dy 

for all X e 0.x, y such that y + h{x) G D^. 

Assumption 4. There exist constants C4, C5, ce, cy, cg, cg, cio, cn, and C12, 
such that 

dVs 

-^Af{x,y + h{x)) < C4||a;||||j/|| 
9Vb , ^^ / „ ,,2 



dx 



■f{x,y + h{x))<C5\\y\\ + cq\\x 



^ A,(a:, y + h{x), e) < e(cr||yf + C8||x||||j/|| + c\\y\\) 

.OVb dh o II II 

+ /i(a;))| < cio||y|| + cii||a;||||y|| + ci2||y|| 

for all X G D^, y such that y + h{x) € Dz- 

These assumptions are equivalent to the assumptions in the classical liter- 
ature, just that here the terms cg||y|| and Ci2||y|| are added. Referring to the 

discussion at the beginning of this section, these terms are relevant for cases 
where g and / are not sector bounded in ||a;|| and ||y||, but may have a bounded 
offset in e. 

First, we show that x will remain within in the perturbed system, so that 
(11) can be used later on. 

Lemma 1. Under assumptions 1-4 above, there exists e' such that, for all 
e < e', the solution x{t) of (1) stays in Cl^ for allt>0 and any initial condition 
where xq S {rf2||a;|p < gD}, with < ^ < 1. 



Proof. From the classical results on singular perturbation [7, Th. 11.2], it can 
be established that the solution of (1) satisfies < 7ie~^ +73£ with some 

positive constants 71,2,3, for t > 0. Considering Vs as a Lyapunov function for 
(1), we find 

= ^ ^^""^^ + y + ^^""^^^ (15) 

< -C2||a;||^ + C4||a;||(7ie~^ +73e)- 

Thus there exists T > 0, such that x{t) e {rf2||a;||^ < D} for t e [0,T] and 

sufficiently small e. For t > T, we have Vs < — f;2|l.'z;|P + C4|la;||74£ and thus 
Vs <0 whenever d2||a;|p = D, for sufficienly small e. This implies that x{t) S 
for all t > 0. □ 

Theorem 1. Suppose that assumptions 1-4 above hold and let 5, 7 such that 

< 7, (5 < 1. Define w(t)'^ = {x{t), z{t) — h{0))'^ , where (x(t), z{t)) is a solution 
of system (1). Let xo G {(i2||a;|p < gD}, with < g < 1. Then there exists 
e*{j,d) > such that for any e < s*{^,5), 

\\wit)f<P{\\w{0)f,t)+Ke), (16) 

where P is a class KC function, i.e., strictly increasing in the first argument, 
strictly decreasing in the second one, and I3{r^, 00) = /3(0,t) = 0, and /x(e) 
as e ^ 0. 

Proof. We first construct a transformed system by letting y = z — h{x), and 
denote the state of the transformed system as = (x,y)^ . From (1), the 
transformed dynamics are 

X = f{x,y + h{x)) 

ey = g{x, y + h{x), e) - ^^fi^, y + h{x)). 
For the transformed system, consider the Lyapunov function candidate 

V{x,y) = il-6)Vsix) + SVB{x,y). (18) 

Let us first check that V satisfies quadratic norm constraints from below and 
above. In fact, from (10) and (12), we directly have 

{l-d)d4x\\^+dd3\\y\\^<V{x,y) 

and 

V{x,y) < {1-S)d2\\x\\^ + Sd4y\\^. 
Next, let us compute the derivative of V{x, y) along trajectories of the system 



(17). Using the bounds in Assumptions 2-4 and Lemma 1, we obtain 
V{x,y) < -(1 - (5)c2||x||^ + (1 - '^)c4||a;||||?/|| +(^0511^11^ 

+ fc6||x||||y|l--C3||?7|p+<5c7||y||2 + <5c8||x||||2;|| 

C 

+ 6cQ\\y\\ +e(5cio||y|p +e(5cii||x||||y|| + edci2\\y\\ 

<-{\\xl\\y\\)Ai^,5,e){\\x\\,\\y\\r 
+ SB{x,y,'y,e) 



with 



and 



(1-5(1-7))C2 _ ^4 + ^(c6+C8-C4+^cii) 

S{c7 + eclo) 



e4 + (5(c6+e8— e4 + eeii) ^£3(1 — y) 
2 



B{x,y,^,s) = (cg +eci2)||y|| - -7C2||a;||^ 

Clearly, if S and 7 satisfy the conditions in Theorem 1, we can find e"{6, 7) > 
such that A{j,S,e) is positive definite for all positive e < £"(7,(5). By con- 
sidering the determinant and trace of A, one can even compute suitable s", but 
the result is a bit longish and skipped here. 

From Lemma 3 in the Appendix, we conclude that B{x, y, 7, e) < whenever 
||u||2 = ||a;||2 + ||y||2 > ^(g-)^ ^j^ij ^ as £ ^ 0. 

Thus, from results on boundedness and ultimate boundedness [2, 7], we can 
conclude that solutions of (17) are bounded by 

Hi)f <^(H0)f,i) + /i(£), (19) 

where /3 is a class ICC function. From Lemma 3, we have /i(e) — >■ as £ — > 0. 

Finally, let us transform the bound (19) back to the original variable = 
{x, z — h{0))'^. From assumption 1 and a Taylor expansion of h{x) — h{0), with 
convexity of D^, we have 

\\wr<\\xr + \\z-hix)r + \\h{x)-h{o)r 

<\\xf+\\yr+ci\\xr 
<{i+ci)\\vr. 

Similarly, 

\\vr<\\xr + \\z-hio)r+\\h{o)-h{x)r 

<(l + c?)lkf- 

Also, note that a <b implies that (3{a,t) < [5{b, t). Thus, we conclude 

\\w{tw<{i+ci){m+ci)\\wm\t)+m)- (20) 

□ 



3 DYNAMICS OF A COUPLED CELL POP- 
ULATION 

3.1 Modeling 

We consider a population of N cells which exchange molecules with a joint 

extracellular medium via diffusion. The concentrations of these molecules within 
cell i are denoted by the vector e . The intracellular dynamics, neglecting 
coupling effects, are assumed to be given by the differential equation 

= F{e% (21) 

The vector field F(^) is assumed to be continously differentiable. 

Next, let us turn to the diffusive coupling. The concentration of the molecules 
in the extracellular medium is denoted by G R^. Diffusive exchange be- 
tween cells and the extracellular medium is described by a flux rate which is 
linear in the concentration difference ^^'^ — (,^^\ with a diffusion rate constant 
kc > 0. Due to volume differences between the intracellular and the extracel- 
lular compartments, a scaling factor needs to be introduced. We will assume 
that all cells have identical volume Vc and that the volume of the medium is 
NVc, i.e., each cell shares a fixed part of the medium equal to its own volume. 
In addition, we include a linear removal rate for the considered molecules from 
the extracellular medium, for example by degradation or an outflow of medium, 
with a rate constant k(i> 0. 

The overall dynamics of the coupled population are then given by the fol- 
lowing system of differential equations. 

^(0) = -fc,^(")+^A^(^O)-^(0)) 
^W=F(^W)_^(^W_^(o)), i = i,...,7V 

3.2 Analysis via singular perturbation theory 

We are going to analyse the population model (22) with the theory developed in 
Section 2. Since we are interested in the limiting case of fast diffusive exchange 
between cells and the extracellular medium, a reasonable small parameter is 
e = k~^. First, we have to bring the model into the standard form given in (1). 
To this end, consider the variable 

^ = ^^°' + ^E^^'^ (23) 

which from (22) is subject to the dynamics 

1 ^ 

i = -ka^«^) + ±Y.n^^'^)- 



(22) 



In order to have the notation as in Section 2, we also denote z^*^ = 
1, . . . , A*". Thus, the transformed model with variables x and z is 

N N 

X = -k,{x '^'^) + ^(^^'^) 



(24) 



which is in the standard form for singular perturbation analysis with small 
parameter e = k~^. 

To simplify notation, we introduce / and g such that (24) writes as 

x = f(x,z) , ^ 

. ) \ 25 
ez = g[x,z,e). 

Next, we derive the slow and boundary-layer models for the population 
model (24) according to (6) and (7). The first step is to solve the equation 
g{x,z,0) = to obtain the quasi-steady state solution z = h{x). From the 
model (24), this yields the system of nN linear algebraic equations 

1 " 

z^"^ -x+-^z^^^ =0, i = l,...,N, (26) 

i=i 

which has to be solved for z^*' in order to obtain h. One can verify directly that 
one solution of (26) is 

^W = |, i = l,...,N. (27) 
To show that (27) is the unique solution, we have the following result. 

Lemma 2. The matrix E^^ = ^1 + / G M^^^, where 1 is a N x N matrix of 

all ones, has one eigenvalue at 2 and N — 1 eigenvalues at 1. 

Proof. First, let us compute the eigenvalues A of ^1. From 

where 1 G is the vector of all ones, we can verify that one eigenvalue is 
Ai = 1. Since ^1 is a rank one matrix, the other eigenvalues are A2...JV = 0. 
The eigenvalues A of Ejsf are computed by 

det((A-l)7-ll) = 0, 

yielding 

Ai = Ai + 1 = 2 
A2...JV = A2...Ar + 1 = 1. 

□ 



Note that (26) is actually a system of n decoupled linear equations in the 

(i) 

N variables zj. , i = 1,...,N, for each of fc = 1, . . . , n. With some abuse of 
notation, we have for each component of z^^^ 



EnzI*'' = Xk- 

Applying Lemma 2, since En does not have zero eigenvalues, we conclude that 
this equation has a unique solution for each component z^. , and thus the alge- 
braic equation (26) has the unique solution given by (27). 

The slow model is then obtained by substituting (27) into f{x,z), yielding 

x = f{x,hix))=F{^)-k,^. (28) 

The boundary- layer model describes the dynamics of the deviation of the 
fast variables z^'^^ from their slow-time-scale approximation (27). As described 
in Section 2, the boundary-layer model is formulated in the variable 

y^'^=z^'^-^. (29) 
Using the model equations (24), the boundary-layer model is obtained as 

'^=-W^Jlty'''^^ ^ = ^,■■■,N. (30) 

Proposition 1. Assume that x is an equilibrium point of the slow system (28), 
i.e., F{x/2) — kixjl = 0. Let the following conditions be satisfied: 



1. F is globally Lipschitz, i.e., ||F(0 - F{C)\\ < L\\^ - C\\ for & 



2. Assumption 2 from above is satisfied for the function f{x, h{x)) = F{x/2)— 

kdx/2 for Dx = K" , where \\x\\ is replaced by \\x — x\\ to account for the 
non-zero equilibrium point, and the Lyapunov function Vs{x) satisfies 



BV 

\^(x)\\<c\\x-x\\ (31) 



for a constant c > 0. 



Define w{t) = i^^^^t) — x/2)i^o,...,N, where ^^*'(i) is a solution of (22). Let 
xo G {d2||a; — < qD}, with < g < 1. Then there exists £* > such that 
for all e < £* , the solutions of the system (22) are bounded by 

\\w{t)\\<l5{\\'^m,t)+n{e), (32) 

and iJ,{e) — )• as e — )• 0. 



Proof. The result is a direct consequence of Theorem 1. In the proof, it remains 
to verify that the transformed system (24) satisfies the assumptions made for 
theorem 1. 

For assumption 1, we have shown previoiisly that the system (24) has a quasi- 
steady state solution z^^^ = x/2, i = 1,. . . ,N. The norm bound in assumption 
1 is thus 

dh 

Prom the assumption in Proposition 1, {x, h{x)) is an equilibrium point of the 
system (24). 

Assumption 2 from Section 2 is a condition in Proposition 1, transformed to 
the non-zero equilibrium point. 

Assumption 3 concerns the boundary layer model (30). This is a linear 
model, and with Lemma 2 we can conclude that its dynamics matrix has n 
eigenvalues at —2/Vc and n{N — 1) eigenvalues at —1/14. It follows that there 
exists a Lyapunov function Vb for the boundary layer model which satisfies 
assumption 3 in Section 2. Also, since the boundary layer model does not 
depend on a;, Vb can be chosen such that OVb/Ox = 0. 

For assumption 4, we make use of the conditions that F is globally Lipschitz 
and thus f,g are globally Lipschitz in x, z, and £, the derivative of Vs{x) is 
bounded by (31), and OVs/dx = 0. Prom these conditions, one can directly 
verify that there exist constants C4, . . . ,Ci2 such that the inequalities (14) are 
satisfied. 

The bound (32) is then a direct consequence of Theorem 1. □ 



4 A POPULATION OF COUPLED SWITCHES 

4.1 Modeling and analysis 

The dynamics of a population of coupled cells become particularly interesting 
for intracellular dynamics with non-linear behavior such as bistability or oscil- 
lations. Here, we focus on bistability, which has been thouroughly studied in 
the context of intracellular biological networks [1, 5, 4]. A typical model for an 
intracellular switch is given by the scalar differential equation [13] 

^ = = if^ " ^' ^^^^ 

where ^ e M denotes the expression level of an auto-activating gene. 

We consider a system where the auto-activating gene product can be ex- 
changed among cells. Using the modeling approach presented in Section 3, the 



resulting population of coupled switches will be described by the system 



^■=1 ' (34) 

3(^W)^ (i) fee (,) (0) 

where the degradation in the extra-cellular medium has been set to zero. 

From the theoretical results in Section 3, one can conclude that the popula- 
tion of coupled switches behaves as a single switch for sufficiently fast diffusive 
exchange: due to the fast coupling, the individual cells would quickly reach a 
state where all have approximatively the same concentration ^^^^ « ^^^^ ss • • • « 
and from then on the dynamics would be governed by the slow bistable 
equation 

x=^^-^. (35) 

1 _|_ 2 ^ ' 

The initial condition of the slow model is computed from the initial condition 
for the full model (33) by the coordinate transformation (23), yielding 

a.(0) = e(°n0) + ^Ee^'H0) (36) 

The inverse transformation 

^ii)^t) = ^, i = 0,...,N (37) 

can be used to obtain an approximate solution of the original population model (34) 
from a solution of the slow model (35). 

Motivated from these theoretical considerations, we perform a simulation 
study with a range of diffusion rate values and heterogeneous initial conditions 
for the switch population (34). Simulation results for A'' = 10 are shown in 
Figures 1 and 3. The solutions for e = have been obtained by simulating the 
reduced model (35) and using (37) to compute the approximate solution for the 
full model. By construction, this approximate solution has identical values for 
all individual switches. 

As indicated by the theoretical considerations, we in fact observe switch syn- 
chronization and population level bistability for the coupled switch model (34): 
for sufficiently small £, the states of individual switches quickly converge to each 
other, and then, depending on the initial conditions, converge all together to 
either the upper or the lower stable equilibrium point. 

In contrast, for larger values of e, individual switches may converge to dif- 
ferent stable equilibrium points, depending on their individual initial condition. 
An example for this is shown in Figure 3 with £ = 20. 
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Figure 1: Simulation of the coupled switch population with iV = 10 for different 
values of e and initial conditions ranging from 0.1 to 1. 
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Figure 2: Simulation of the coupled switch population with iV = 10 for different 
values of e and initial conditions ranging from 0.03 to 0.3. 
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Figure 3: Simulation of the coupled switch population with iV = 10 for different 
values of e and initial conditions ranging from 0.1 to 4. 

4.2 Discussion of biological relevance 

As shortly discussed in the introduction, the relevance of bistable dynamics in 
gene regulation and biochemical signal transduction is very well established. 
However, the relation between bistable intracellular dynamics and population 
level behavior of coupled cells has not been well studied so far. 

Motivated from the theoretical results in Sections 2 and 3, the simula- 
tion study described in this section indicates the occurence of population-level 
bistablity for fast enough diffusive exchange. We speculate that this type of 
coupling may be a relevant biological mechanism to provide robustness of cellu- 
lar decisions within a cell population under cellular heterogeneity and molecular 
noise. 

Coordination among bistable switches in a cell population may for example 
ensure that a cluster of cells within a developing tissue differentiates together, 
even if not all cells receive a high enough stimulus to initiate single cell difer- 
entiation. Alternatively, this mechanism may be useful to extend engineered 
genetic switches on the single-cell level [6] to coordinated switching of a whole 
cell population using synthetic biology techniques. 

5 CONCLUSIONS 

We have proven convergence to a common equilibrium point for populations of 
cells which are coupled via diffusive exchange of molecules with a common ex- 
tracellular medium, under the assumption of sufficiently fast diffusive exchange 
rates compared to the intracellular dynamics. While we have considered only 
homogenous intracellular models in this paper, the results extend directly to 
heterogeneous cells. 

A study of coupled bistable switches showed coordinated behavior and population- 
level bistability among the switches. This has potential relevance for the robust- 



ness of tissue differentiation and for synthetic biology. 



APPENDIX 

Lemma 3. Let 

B{x,y,e) = Ue)\\y\\ - kM'' - (38) 
with k\, k2, e > 0. // 

+||y|| > ^2 + 4fc^fc3 . (-^y) 

then B(x, y, e) < 0. 

Proof. Case 1 Let \\yf > Then 

B{x,y,e)<k,{e)\\y\\{l-Jj^^M)<0. 

Case 2 Let ||y||2 < Then ||a;|p > and, by quadratic exten- 

sion, 

B{x,y,e) = -k4xr - '^{hf - ^-^f + 



^3 

□ 
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